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NUMERICAL SIMULATION OF THE FLOWFIELD 
OVER ICE ACCRETION SHAPES 


1 . INTRODUCTION 

This is the first Semi-Annual Status Report submitted on 
Grant NAG 3-665. It includes the progress made during the period 
of October 1985 to April 1986. The NASA Technical Officer for this 
grant i 3 Dr. Robert J. Shaw, NASA Lewis Research Center, Cleveland, 
Ohio. 

2. R&D STATUS REPORT 

The grant program schedule consists of five components, i.e., 
computer code revision, grid generation, computations, extraction 
of data, and report preparation. During thi3 reporting period the 
first two components were accomplished and progress made on the 
computations. In addition, an abstract was prepared and submitted 
to the AIAA for potential presentation of these results at the 
Aerospace Sciences Meeting, January 12-15, 1987. Details of the 
program status follow in the next sections. 

2 . 1 Technical Effort 

In this section the progress of the technical effort 
will be discussed. The primary goals of this program are directed 
toward the development of a numerical method for computing flow 
about ice accretion shapes and determining the influence of these 
shapes on flow degradation. In pursuing these goals, it is 
expedient to investigate various aspects of icing independently in 
order to assess their contribution to the overall icing phenomena. 

The specific aspects to be examined include the water 
droplet trajectories with collection efficiencies and phase change 
on the surface, the flowfield about specified shapes including 
lift, drag, and heat transfer distribution, and surface roughness 
effects. In treating these issues it must be kept in mind that 
they will ultimately be coupled together in a fashion which will 
permit accurate prediction of the complete icing process as well as 
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the influence on the surrounding flowfield. The following 
paragraphs describe the progress in examining these issues. 

The first issue to address is the unsteady nature of the 
phenomenon . 1 Although the ice shape configuration changes with 

-I4 

time, the rate of growth is slow (10 ft/sec) compared to the 
flight speed (200 ft/sec), therefore, a quasi-steady analysis is 
appropriate. This approach involves computation of Navi er-Stokes 
solutions for a sequence of shapes starting from the original 
configuration, computing the growth rate, advancing the 
configuration over a fixed time interval (e.g., 50 sec), 
recomputing a new growth rate, etc., until the final period is 
attained. 

A primary key to this approach is the accurate computa- 
tion of the heat transfer over a roughened surface of unusual 

geometry. Although the Navier-Stokes equations conceptually 
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possess the ability to simulate this process, numerical 

confirmation has not yet been achieved. Fortunately, a series of 

5 i) 5 6 

experimental heat transfer tests have been accomplished ’ ’ 
which assist in the validation of the computational fluid dynamic 
(CFD) methods. The purpose of this phase of the investigation is 
the computation of a series of these test cases to validate the 
numerical simulation. 

The configurations computed were models of ice accretion 
shapes formed on a circular cylinder in the NASA Lewis Icing 
Research Tunnel (IRT). These shapes were 2-, 5-, and 15-minute 
models of glaze ice and a 15-minute accumulation of rime ice. An 
existing Navier-Stokes program was modified to compute the 
flowfield over these four shapes. 

2.1.1 Governing equations 

The governing equations are obtained by adding 
body forces to the Navier-Stokes equations to account for the drag 
of the droplet particles and surface roughness. In addition, a 
continuity and momentum equation are required to develop the 
trajectory equations for the droplets. 
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where 

I ■ -55" IS • Spl<Y - V p ) - droplet drag/unit mass 

The water droplets are assumed to be uniform spheres of diameter, 
d. Cloud physics experiments indicate d to vary between 5 and 50 
microns depending upon the air temperature and liquid water content 
(see Aircraft Icing, NASA CP2086). Often an average value of 20 
microns is used in ice accretion studies. Therefore, 
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C D = C D (Re d ) (see Figure 1) 
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A curve fit of this classic plot produces the 

following: 

Re d < 1 C D 

1 < Re . < 400 C n 

d D 

400 < Re d < 3 x 10 5 C D 

These values will be used in the computation of the water droplet 
trajectories . 

2.1.2 Boundary conditions 

Far field . At the far field boundary, the 
following conditions are prescribed: 

V , T , P 
- 00 00 J 00 

Pp = LWC = liquid water content of cloud 

V = V 

- P -00 

Surf ace . At the surface, the following 
conditions prevail: 

V = 0, T specified = 0 

— w 9n 

2.1.3 Grid generation 

The grid was developed using the hyperbolic grid 
generator of Steger and Barth. This technique appears ideally 
suited for computing grids for the irregular shapes of icing 
configurations. This method first requires a detailed description 
of the surface geometry. The program produces an orthogonal, body 
oriented grid (which is preferable in CFD computations) and 
clusters the grid points near the surface. The hyperbolic method 
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(Stokes flow) 
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= 0.5 
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has little control over the final 
perfectly acceptable for external 
for internal flows). An example 


exterior grid shape, which is 
flow problems (but unacceptable 
grid is shown in Figure 2. 


2. 1.3*1 Solution of the water droplet 
equations 


An efficient means for solving the 
water droplet trajectory equations was developed during this 
reporting period. Droplet trajectory equations in the past have 
been solved using the Lagrangian method. In this investigation it 
was proposed to use the Eulerian method for two reasons. First to 
make the droplet system of equations compatible with the airflow 
equations, and secondly to include the variation of the liquid 
water content (P p ) throughout the flowfield. The Lagrangian method 
presumes the water content to be constant. 


The droplet equations are similar in 
form to the airflow equations with the exception that the stress 
tensor is zero. This difference, however, changes the mathematical 
character of the partial differential equations from elliptic to 
hyperbolic as will be shown in the following paragraph. This 
observation means that the hyperbolic equations may be marched in 
space using a Parabolized Navi er-Stokes code (PNS). An advantage 
results in that the droplet equations may be solved in only a few 
seconds on a VAX computer. 

2. 1.3. 2 Eigenvalue analysis of marching 
techni que 


The purpose of this section is to 

demonstrate the applicability of space-marching techniques to solve 
the governing fluid dynamic equations for water droplets motion. 

To correctly apply a space-marching technique the eigenvalues of 
the governing equations must be real (indicating the character of 
the equations is hyperbolic). 


form are written as: 


The governing equations in divergence 


E + F = H 
x y 
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Figure 2. Grid for 15-Minute Glaze Case. 
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where E, F, and H are defined by 
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In order to perform the eigenvalue 
analysis the E and F vector are factored such that 

E = § U = A U 
x 3U -x - -x 

£ F 

F = ~ U = B U 
y 3U -y - -y 


where 
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The A and B matrices are found to be 
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Applying the above factorization can be rewritten as 


U + A -1 B U = A _1 H 
-x - - -y 
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To determine the character of this equation set it is necessary to 
find the eigenvalues of the matrix A ^B. A ^ B was found to be 


A _ 1 B = 


v 
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0 


2 

u 


v 

u 


p_ 

u 


0 

V 

u 


-1 


The characteristic equation of the 


matrix A B is determined from 


det[A _1 B - XI] = 0 


or in simplified form 




= 0 


The eigenvalues of the matrix A 1 B are defined as the roots of the 
characteristic equation. Therefore, the eigenvalues of A 1 B are 


X 


_v 

u 


Since the eigenvalues are real, the character of the governing 
fluid dynamic equations is hyperbolic. Therefore, the governing 
equations can be solved with a hyperbolic space-marching scheme. 


2.1.4 Solution procedure 


The airflow was computed using a well documented 
1 4 

Navier-Stokes code adapted from the original program developed by 
Shang. The MacCormack algorithm is utilized and the program is 
vectorized for the CRAY computer. Also, the hyperbolic grid 
generator of Steger and Barth was used as discussed in the previous 
section. Since this portion of the project is current state of the 
art in CFD , little difficulty was encountered. 


The computation of the droplet trajectories is 
computed using technology developed for the. PNS solving schemes. 
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The droplet equations are hyperbolic since the stress- tensor 
vanishes for this model. (No collisions between droplets are 
considered in the formulation and hence no pressure or shear terms 
are present.) The governing equations are transformed into a £-n 
computational domain and take on the following form: 
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where the contravari ant velocity components are represented as 


U = y u - x v 
p n P 
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The upstream boundary conditions are 


p p («) 

y-> 

Vp (.) 


LWC 

V 

00 

0 


Given the air velocity components (u,v) from the solution of the 
airflow equations, the droplet equations are marched in the n~ 
direction inward towards the body. After attaining the body 
surface the values of Up,Vp, and Pp are recorded from which the 
local collection efficiency (B) may be determined. 
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where tan 8 


/W \ 

ydx/w 


= slope of the wall. 


A dimensional analysis of the water droplet 
equations indicates that the collection efficiency is a function of 
the geometry and only one flow parameter, t’; where 


This single parameter is a type of Reynolds number, taking into 
account all flow variables of the problem, and is very useful for 
evaluating collection efficiencies. Figure 3 shows the sensitivity 
of 8 to t ' for various flows over circular cylinder. 


2 . 2 Results 

To calibrate the numerical procedure the heat transfer 
distribution over a smooth circular cylinder under laminar condi- 
tions was computed for which there exists an exact solution 
•• 1 5 

(Frossling ). Figure 4 shows a comparison of Nusselt number over 
a cylinder for the Navier-Stokes solution with the Frossling 
solution. Conditions for this case are Re^ = 138,000 and a 
diameter of 2 inches. Excellent agreement may be observed, proving 
the validity of the numerics. 

Figure 5 shows the computed droplet velocity vector field 
around a 2-inch diameter cylinder in potential flow with a free 
stream velocity of 130 fps. 

Preliminary computations of the flowfield about each of 
the ice shapes have been performed. Refinement of these computa- 
tions are currently being examined in order to achieve the required 
accuracy . 
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COMPUTED COLLECTIVE EFFICIENCIES 
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Figure 3. Computed Collection Efficiences for Various x' Values for Cylinder. 
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Figure Comparison of Computed Nusselt Number 

Distribution with Frossling Solution. 
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COMPUTED DROPLET VELOCITY VECTORS 



x 

Figure 5. Computed Droplet Velocity Vectors 
Over a Cylinder x' = 0.476. 
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APPENDIX 


Description and Use of HGRID 


Hyperbolic Grid Generator 


HGRID has a nicely structured modular form. Though it is not 
completely clear to the uninformed observer ju3t what each subroutine does, 
the important routines are easily identifiable. Much of the code should be 
considered a "black box". The routines which you will most likely want to 
change are the main program, BODIS, and INITIA. These sections of the 
program involve setting up the body, defining important parameters such as 
loop sizes and output of the grid. A list of the primary subroutines and 
their purpose is shown below in order of appearance in the code. 


MAIN 


Output of the completed grid. 


SARC - Called at each arc by STEP, sets up integration parameters 

based on input parameters and scaling factors. 

BODIS - Reads in or generates the body coordinates, called by 

INITIA. 


INITIA 


METRIC 

STEP 

VOLUME 


Reads in or assigns the smoothing and integration input 
parameters, sets up stretching and scaling factors based 
on the input parameters, called by MAIN. 

Generates the coordinate transformation metrics for the 
integration, called at each arc by STEP. 

Marches the grid generation out from the body. Calculates 
X, Y coordinates at next arc, called by MAIN. 

Calculates the "volume" or area of each cell formed by the 
present arc and where the new arc should lie based on 
variable scaling factors. Called by STEP. 


The basic procedure in the process of grid generation is as follows: 

1) Body coordinates and input parameters are read, and the 
body coordinates are written to the output file. 

2) STEP is called and it calls routines to set up: 
integration parameters, cell areas and the right hand side 
of the PDE to be solved. 

3) The PDE is solved by routines called by STEP, and within 
STEP the X, Y values for the next arc are calculated. 
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4) In MAIN, the X, Y values of the new arc are written to 
tape, and STEP is called again. 

The sequence of calculating the next arc out from the body and then 
writing it to tape eliminates the need to store X and Y in two dimensional 
arrays. This cuts down on the memory required by HGRID. 

Changes to HGRID will most likely be changes in the body coordinates, 
loop sizes and final grid form. These changes will all take place in MAIN, 
BODIS and INITIA. The other input parameters defined in INITIA should be 
considered constants that need not be changed. They have been found to work 
well for many grid generation cases. Two factors, however, that you may wish 
to control are DSETA and ESCAL. DSETA is the initial step size out from the 
body in eta, the transformed coordinate normal to the body. ESCAL controls 
the scale factor that is used to calculate the cell areas or volumes. 

Hint: In the event that only the upper or lower half of a grid is needed 

due to symmetry conditions, you must supply HGRID with the whole 3et of body 
coordinates anyway. For HGRID to function, it needs to be given a complete, 
closed body. You can cut the mesh in half after the complete grid has been 
generated. 
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